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Metallic nano-structures characterised by multiple geometric length scales support low-frequency 
surface-plasmon modes, which enable strong light localization and field enhancement. We suggest 
studying such configurations using singular perturbation methods, and demonstrate the efficacy 
of this approach by considering, in the quasi-static limit, a pair of nearly touching metallic nano¬ 
spheres subjected to an incident electromagnetic wave polarized with the electric field along the line 
of sphere centers. Rather than attempting an exact analytical solution, we construct the pertinent 
(longitudinal) eigen-nrodes by matching relatively simple asymptotic expansions valid in overlapping 
spatial domains. We thereby arrive at an effective boundary eigenvalue problem in a half-space 
representing the metal region in the vicinity of the gap. Coupling with the gap field gives rise to a 
mixed-type boundary condition with varying coefficients, whereas coupling with the particle-scale 
field enters through an integral eigenvalue selection rule involving the electrostatic capacitance of 
the configuration. By solving the reduced problem we obtain accurate closed-form expressions for 
the resonance values of the metal dielectric function. Furthermore, together with an energy-like 
integral relation, the latter eigen-solutions yield also closed-form approximations for the induced- 
dipole moment and gap-field enhancement under resonance. We demonstrate agreement between the 
asymptotic formulas and a semi-numerical computation. The analysis, underpinned by asymptotic 
scaling arguments, elucidates how metal polarization together with geometrical confinement enables 
a strong plasmon-frequency redshift and amplified near-held at resonance. 


I. INTRODUCTION 


for the Drude dielectric function is Ilf; 


Advances in micro-fabrication technology along with 
new innovative ideas in nano-photonics and meta¬ 
materials have rekindled interest in the optics of small 
metal particles and structures M- Nano-metallic con¬ 
figurations exhibit unique optical characteristics, which 
are largely associated with the excitation of localized 
surface-plasmon (LSP) modes: conjoint resonant oscilla¬ 
tions of electron-density and electromagnetic fields. LSP 
resonances are manifested in spectral peaks of the optical 
cross-sections and in significant near-field enhancements. 
Such phenomena are attractive for sensing applications 
8 , plasmonic manipulation of light below the diffraction 
limit [5J, enhanced Raman scattering [6j, [7|, plasmonic 
rulers lasers [9f and nanoantennas flol-[l2j. second 
HU and higher-order El harmonic generation, and pho- 
tovoltaics [la ]. 

LSP resonance is exemplified by light scattering from 
a nano-metric metal sphere. In the quasi-static limit, 
the near-field, which appears to be driven by a uniform 
time-harmonic field, is derived from an electric potential 
governed by Laplace’s equation. The induced electric- 
dipole moment is readily found to be proportional to the 
Clausius-Mossotti factor (e — l)/(e + 2), where e is the 
frequency-dependent dielectric function of the metal par¬ 
ticle relative to its dielectric surrounding. The evident di¬ 
vergence at the Frohlich value e = — 2 is associated with 
the existence of a solution to the homogeneous quasi¬ 
static problem (in which the incident field is absent). At 
visible and ultraviolet frequencies, a good approximation 


Re[e] w 1 - u 2 /u 2 , Im[e] « u^y/w 3 , 


(1) 


where ui is the frequency of light, and ui p and 7 are 
the plasma and collision frequencies, respectively (rep¬ 
resentative values for noble metals: Hlo p ~ 9eV, fiy ~ 
0.13 eV). Owing to the smallness of y/w p , a damped res¬ 
onance occurs when Re[e] « — 2. 

LSP frequencies and the magnitude of resonant en¬ 
hancement can be “designed” by shape and material 
variations E3N13- Appreciable modification is achiev¬ 
able with “singular” configurations characterised by non¬ 
smoothness or multiple-scale geo metries: sharp corners 
[U [22j]; elongated particle s [23^ 251: tips (26j; nearly 
touching [27l [20t , touching [3 g T 311, and embedded [32f 
particle pairs; “plasmonic molecules” [13; and particle- 
wall configurations [13 ■ A recurring feature encountered 
in multiple-scale geometries is the existence of modes, 
that with increasing scale disparity, significantly redshift 
to low frequencies and large — Re[e] values. 

Singular configurations of the latter sort, which enable 
a particularly strong near-field enhancement [34|, have 
been extensively studied through numerical simulations 
[35H4(J - analytic solutions in separable coordinate sys¬ 
tems |2ll. l30ll4l[ , and the powerful method of transforma¬ 
tion optics (TO) f42l. [43il . In two dimensions TO furnishes 
explicit solutions 22 T 2 ll [28j that can be generalised to 


address effects of retardation [44] and non-locality [45[ . 
TO is less straightforward to apply in three dimensions, 
and typically furnishes an efficient computational scheme 
rather than simple analytic expressions; neverthe¬ 
less, the method has been successfully applied in vari¬ 
ous scenarios, e.g. in the design of broadband plasmonic 
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structures [31] and in the study of van der Waals forces 
[ 20 t . Yet another approach, motivated by general scaling 
arguments [53], is based on approximate reductions of 
implicit infinite-series solutions for nano-wire and nano¬ 
sphere dimers [48l450j ; we delay a critical discussion of the 
accuracy of the expressions in Refs. El and [53 to 9VII 

The above methodologies are “top-bottom” — neces¬ 
sarily relying on “exact” analytical or numerical solu¬ 
tions, which then need to be reduced or evaluated in sin¬ 
gular geometric limits (e.g. small gaps for dimers, increas¬ 
ing aspect ratio for elongated particles). In this paper, 
we suggest an alternative “bottom-up” approach, which 
is specifically suitable for analyzing the extreme phys¬ 
ical characteristics of LSP resonance in multiple-scale 
metallic nano-structures. The approach is built upon 
the paradigm of singular perturbation methods, where 
an asymptotic rather than an exact solution is sought by 
carrying out an ab initio singular analysis. In particu¬ 
lar, applying asymptotic scaling arguments in conjunc¬ 
tion with the method of matched asymptotic expansions 
[5TI ] allows exploiting the spatial non-uniformity of the 
singular LSP modes towards systematically decompos¬ 
ing the problem governing these into several consider¬ 
ably simpler ones. The main outcomes are asymptotic 
expressions for the resonance frequencies, along with a 
minimalistic description of the corresponding modes. In 
the context of excitation problems, when losses are low 
an asymptotic approximation for the field enhancement 
may be obtained with little additional effort. 

To demonstrate this approach, we shall focus in this 
paper on the prototypical configuration of a pair of nearly 
touching metal nano-spheres [H, E3, IH, [53, El, IH, Elr 
[57[| • We limit ourselves to a quasi-static treatment of 
the near field, and focus on the axisymmetric and lon¬ 
gitudinal (bonding) modes that singularly redshift with 
vanishing gap width, and the excitation of these by a 
long-wave electromagnetic plane wave polarized with the 
electric field along the line of centers. The paper proceeds 
as follows. In [JIT] we formulate the problem and note use¬ 
ful generalities. In fJTTT]we consider the plasmonic eigen¬ 
value problem and obtain closed-form expressions for the 
dielectric-function resonance values. In : 1I VI we employ 
the eigenvalues and eigen-potentials found in [JTTT] to ob¬ 
tain formulas for the induced dipole moment and field 
enhancement in the gap under resonance conditions. In 
m we present a validation of our results against a semi- 
numerical solution of the quasi-static problem. Lastly, 
in 1VII we recapitulate our results, compare with related 
works, and discuss recent developments along with future 
directions. 


II. PRELIMINARIES 
A. Problem formulation 
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FIG. 1. Schematics of the quasi-static excitation problem. 


by their radius a, and their frequency-dependent dielec¬ 
tric function relative to the background material, ege, 
6 q being the vacuum permittivity. The minimum sep¬ 
aration between the spheres is denoted by 2 ha. The 
pair is subjected to an incident time-harmonic electro¬ 
magnetic plane wave of frequency w, polarized such that 
the electric field is in a direction i parallel to the line of 
sphere centers, see Fig. [T] We shall focus on the quasi¬ 
static limit of sufficiently small particles, a -C c/ (wy/\e\), 
c being the speed of light in the surrounding dielectric 
[la . l58| . In this description, the pair experiences the inci¬ 
dent wave as a uniform time harmonic field iEYie[e~' lu,t }, 
where E is a real constant magnitude, and the electric 
near-field is approximately curl-free, derived from a time- 
harmonic potential aERe^e - *^]. 

Henceforth, we normalise all lengths by a. In partic¬ 
ular, we denote by x the dimensionless position vector 
relative to the centre of the gap. It is convenient to em¬ 
ploy a cylindrical coordinate system (r, z, 4 >), with origin 
at x = 0 and a 2 -axis directed along i (see Fig. [2]). The 
problem at hand is to find a continuous potential p that 
satisfies 


V • (£Vp) = 0, (2) 

where £ is unity in the dielectric medium and e in the 
metal spheres, together with the far-field condition 

^~-z + o( 1) as |x| —»■ 00 . (3) 

The additive freedom of tp is removed in j3j by choosing 
the disturbance from the incident field to attenuate at 
large distances. 


B. General properties 

The following generalities can be inferred from @ 
([3]), see also Refs. [1, [3f| 15914611 ] . First, the potential 
is axisymmetric about the i axis, i.e. dp/dcj) = 0. In 
addition, the potential is antisymmetric about the plane 
z = 0. It is therefore possible to consider only the half¬ 
space z > 0 , with 


Consider an identical pair of metal spheres surrounded 
by a dielectric medium. The spheres are characterised 


p = 0 at z = 0. 


( 4 ) 
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We next note that each of the spheres is overall charge 
free. That is, 

/S" = tl - (5) 

where the integration is over an arbitrary closed surface 
entirely within the dielectric medium, which encloses one 
or both spheres. Thus p lacks a monopole term in its far- 
held expansion 0 , which can therefore be extended to 

£ 

p ~ — 2 + /i-j— pj as |x| —> oo. ( 6 ) 

l x r 

The complex-valued constant p is the induced dipole mo¬ 
ment of the pair, normalised by 4 - 7 ra 3 eo edE (e,i being the 
relative dielectric constant of the surrounding medium). 

Green identities applied to p and its complex conjugate 
p* , in conjunction with 0 and ([ 6 ]). can be shown to give 

47 rlm[//] = Im[e] f V<p • Vp* dV, (7) 

J metal 

where the integration is over the volume of the two 
spheres. Eq. 0 is a work-dissipation relation: The work 
done by the incident held on the effective induced dipole 
of the dimer is lost to Ohmic dissipation in the metal. 
This relation serves to remind us that we are in the quasi¬ 
static limit, where radiation losses are negligible. We 
note that 4 - 7 rIm[/x] gives the quasi-static approximation 
for the absorption cross-section normalised by toa 3 /c 
Excitation problems of the above sort are ill posed for 
a set of discrete, real and negative, e values, for which the 
homogeneous variant of the problem obtained by omit¬ 
ting the incident held in 0 possesses non-trivial solu¬ 
tions. Since in practice Im[e] ^ 0, such solutions are 
often termed “virtual” LSP modes. Consider now an 
eigen-potential p r corresponding to one such eigenvalue 
e r . At large distances, 

p r ~ p r — 13 as |x|->■ oo, ( 8 ) 

l x l 

where p r is arbitrary. We note that the integral relation 
0 does not apply to p r , since condition 0 is used to 
derive it. Say, however, that p is the solution to the 
actual quasi-static problem, with e = e r + Ae. Applying 
the Green identities to the pair of potentials (p, p r ), one 
can easily show that 

4:np r = — Ae / Vp r • VpdV. (9) 

J metal 

Whereas 0 bears a clear physical signihcant, we shall 
find in IIVI that 0 is more useful for deriving quantita¬ 
tive estimates for the resonant held enhancement. 


C. Near-contact limit 

Henceforth we consider the near-contact limit where 
h —► 0. Specifically, our goal will be to derive closed-form 
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FIG. 2. Schematics of the dimensionless excitation problem. 


asymptotic formulas for the eigenvalues e^ n \ and the val¬ 
ues at resonance conditions (i.e. when Re[e] = e("') of the 
induced-dipole moment p and the held enhancement in 
the gap, G = - (dp/dz) x ^ Q . 

III. LSP MODES 

In this section we consider the boundary eigenvalue 
problem obtained by replacing the far held 0 with the 
homogeneous condition p —> 0 as |x| —> oo. It is clear 
that the only relevant modes are those that are both 
longitudinal, i.e. anti-symmetric about the plane z = 0 , 
and axisymmetric. 


A. Scaling arguments 

How do the eigenvalues e scale with h7 In the vicinity 
of the gap the spherical boundaries can be approximated 
by two oppositely facing paraboloids. The separation 
therefore remains 0(h) for 0(h 1 ^ 2 ) radial distances. We 
shall refer to the subdomain of the dielectric medium 
where z = 0(h) and r = 0(h 1 ^ 2 ) as the gap domain, 
and assume without loss of generality that the eigen- 
potential p corresponding to e is 0 ( 1 ) there; since p is 
antisymmetric about the z plane, the transverse gap field 
must be 0(1/h). Consider now the particle-scale dielec¬ 
tric “outer” domain, away from the gap, where p varies 
over distances |x| = 0(1) and is at most 0(1). It follows 
from 0 that p = 0(h) when 2 = 0(h) and r = 0 ( 1 ). 
This implies that the gap potential must vary in the ra¬ 
dial direction over O ( h x / 2 ) distances. Consider next the 
metal sphere in z > 0 , in the vicinity of the gap boundary. 
Continuity of potential implies that p is 0(1) there, and 
varies rapidly over O ( h : / 2 ) distances along the bound¬ 
ary. On the latter length scale, the metallic volume ap¬ 
pears transversely semi-infinite. Thus, the symmetry of 
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Laplace’s equation implies comparably rapid transverse 
variations. We shall refer to the metal subdomain where 
r,z = 0{h 1 ^ 2 ) as the pole domain. 

The scaling of e now follows by considering the continu¬ 
ity of electric displacement across the gap-pole interface: 
The product of e and the 0(h -1 / 2 ) transverse field in the 
metal pole must be on the order of the 0(1 /h) transverse 
field in the gap. We accordingly find 

e~-h~ 1/2 \(h) + 0(l), 0 < A = 0(1), (10) 

where, as we shall see, the pre-factor A itself depends 
weakly upon h , being a function of ln(l /h). This more 
subtle feature can also be anticipated from scaling argu¬ 
ments, as we now discuss. 

The largeness of |e| suggests that, on the particle- 
scale, the spheres are approximately equipotential. Thus, 
the particle-scale potential distribution in the dielectric 
medium is approximately that about a pair of nearly 
touching spheres held at uniform potentials, say ±V. 
But the latter distribution implies a net charge on each 
sphere, in contradiction with (|5|. This means that the 
deviation in the vicinity of the gap of <p from the particle- 
scale distribution must be consistent with a localized ac¬ 
cumulation of charge such that the overall charge of each 
sphere is zero. In fact, the contributions to the integral in 
© from the gap and outer dielectric domains are readily 
seen to be comparable: The former is given by the inte¬ 
gral of an 0(l//i) transverse field over an 0(h) surface 
area, while the latter is obviously 0(1). This implies a 
connection between the solution in the gap and pole “in¬ 
ner” domains, and the capacitance of a pair of perfectly 
conducting spheres held at fixed potentials ±V, which is 
known from electrostatic theory to depend logarithmi¬ 
cally upon h as h —>• 0 js§|- 

B. Gap and pole domains 



FIG. 3. The resonant modes are constructed in TOIl bv match¬ 
ing simple asymptotic expansions describing the gap, pole, 
and particle-scale “outer” domains. 

Integration of (fl4l) in conjunction with (© gives 

<1- = A{R)Z, (15) 

where the function A(R) remains to be determined. 

Consider next the pole domain of the metal sphere 
in z > 0, described by the alternative pair of stretched 
coordinates 


We begin constructing the eigen-potentials by consid¬ 
ering the gap domain. To this end, we introduce the 
stretched coordinates 

R = r/h 1/2 , Z = z/h. (11) 

The gap domain is where R,Z~ 0(1), with R > 0 and 
Z bounded by the spherical surfaces (see Fig. E| 

Z ~ ±H(R) + 0(h), where H(R) = 1 + ^R 2 . (12) 

As in i jlll A[ we assume without loss of generality that 
<p = 0(1) in the gap domain. We accordingly pose the 
asymptotic expansion 

<p~$(l?,Z) + o(V /2 ) . (13) 

Substituting m into Laplace’s equation, written in 
terms of R and Z, yields at leading order 


R = r/h, Z = z/h 1 / 2 . (16) 

The pole domain is where R, Z ~ 0(1), with R > 0 and Z 
bounded from below by the sphere surface Z ~ Oft 1 / 2 ); 
under the rescaling (fl6l) . the pole domain is to leading 
order a semi-infinite half space. We expand the pole- 
domain potential as 


<p ~ $(R, Z) + O (V /2 ) , 


where $ is governed by Laplace’s equation, 


RdR\ dRJ 



(17) 

(18) 


Expanding the eigenvalue e according to m, requir¬ 
ing leading-order continuity of potential and electric dis¬ 
placement yields 


< 9 2 $ 

5Z2 


= 0. 


9 $ 9 $ 

< 14 > * = *■ az^az ‘ 0 


at Z = H(R), Z = 0. (19) 
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Eqs. © may be combined to yield the mixed-type 
boundary condition 


<9$ 

$ + A H-= =0 at Z = 0, 
aZ 


( 20 ) 


which involves $ alone. We also find the relation 

MR) = jj(*)z = o, (21) 


For p we require a solution to Laplace’s equation that 
attenuates at large distances and satisfies p = V at £ = 1 
and p = 0 at £ = 0. In addition, the solution must 
asymptotically match the gap and pole potentials. Stan¬ 
dard Hankel-Transform methods yield the requisite po¬ 
tential as [111 

p ~ V(£ 2 + rj 2 ) 1/2 [ e -s smh ( 5 ^ J 0 (sri) ds, (26) 
J o smh s 


which provides the gap potential once the pole poten¬ 
tial has been determined. As discussed in ^ 1111A1 the 
particle-scale metal domain is equipotential to leading or¬ 
der. Asymptotic matching thus furnishes the “far-field” 
condition 


-+ V as R 2 + Z 2 oo, (22) 

where V is the above-mentioned uniform potential of the 
sphere in 2 > 0. For later reference, we note that (I2T1) 
and (1221) imply <f> ~ 2 VZ/R 2 as R —> oo. 

The pole potential $ is governed by Laplace’s equation 
m in the half space Z > 0, the boundary condition 
and the far-field condition (1221) . We shall see in 
^III FI that this boundary-value problem possesses non¬ 
trivial solutions for essentially arbitrary A > 0 and V, 
excluding a discrete set of A values for which a solution 
exists only if V is set to zero. It is tempting to guess 
that the latter set furnishes the requisite eigenvalues; we 
shall see in fllll El however that these only give a leading 
term in an expansion of A in inverse powers of ln(l /h). 
Avoiding such gross logarithmic errors, the eigenvalues 
cannot be solely determined from local considerations, 
as already anticipated in i llll Al We accordingly divert at 
this stage to consider the particle-scale dielectric domain. 


C. Outer domains and global constraint 

In the outer limit, h -+ 0 with x = 0(1), the spheres 
appear to be in contact to leading order. It is therefore 
convenient to introduce the tangent-sphere coordinates 
defined through [63} 


r+ »r 


(23) 


The plane z = 0 is mapped to £ = 0, the upper-half-space 
sphere to £ ~ 1 + 0(h), and the origin to 77 —► 00 (see 
Fig. I3J). We expand the potential in the outer dielectric 
domain as 


p-p + Oih 1 / 2 ), (24) 

and in the outer metal domain as 

p~ ±V + 0(h 1/2 ), (25) 

the ± sign corresponding to the sphere in z sj 0, respec¬ 
tively. 


where Jo is the zeroth-order Bessel function of the first 
kind. It may be readily shown that p ~ 2 Vz/r 2 as 
77 — > 00, whereby (1261) matches with the gap and pole 
domains. Asymptotic evaluation of (l26l) for small £ and 
77, together with C51) . yields the asymptotic behaviour of 
(l26l) at large distances [cf. ©], 

, 7T 2 V 2 

<P - as l x l —^ °°- (2 ? ) 

The leading-order matching of (l26l) with the solution 
in the gap and pole regions has been achieved, appar¬ 
ently without introducing any additional constraint on 
the problem governing $. Instead of continuing the 
asymptotic analysis to the next order in search of the 
awaited eigenvalue selection rule, it proves simpler to 
follow the discussion in £1111 Al and consider the charge 
constraint dSJ) applied to the sphere in z > 0. In accor¬ 
dance with the discussion there, an attempt to calculate 
the contributions to the integral in ([3|) from the gap and 
outer region separately results in diverging integrals 
this is typical in cases where contributions from overlap¬ 
ping asymptotic domains are comparable (5l| . Following 
the analysis in Ref. [(32}, we write (O as 



A(R)RdR + 



dp 


T^ d ’ 1=0 { hI ' 2 )’ 

(28) 


where 770, lio> 1 are connected through = 2/770+ 

0(1/?7q). The first integral in (EBl) . which represents the 
excess contribution from the gap region, can be written 
as 



V 

WV_ 


RdR + V\nH(R 0 ). 


(29) 


The second integral in (I28[) . which represents the contri¬ 
bution of the outer domain, can be worked out to be 


~ V [2 7e + ln(l + 7T§)] , (30) 


where y e ~ 0.5772 is the Euler-Gamma constant. Noting 
that \nH(Ro) ~ ln(2 /h) — ln77g, and invoking fl!HT) . the 
singular parts of the two contributions cancel out, and 
we find the anticipated “global” constraint: 


ln(2 /h) + 2y e 



$(R,0) - V 
H(R) 


RdR. 


(31) 
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D. Reduced eigenvalue problem 


Together with (1311) . the problem governing the pole 
potential $ constitutes a reduced boundary-eigenvalue 
problem. It is convenient to recapitulate this problem in 
terms of the disturbance potential, 2) = <f> — V. We 
have Laplace’s equation 


RdR\ dRj 



(32) 


in the upper half-space Z > 0, the integro-differential 
varying-coefficient boundary condition 


0 + Atf(i?)|| 


1 

ln(2 /h) + 2q e 



fjR, 0) 

H(R) 


RdR (33) 


where Vo(s) depends on the conditions at 2 = 0. Trans¬ 
forming the boundary condition (13611 , we find the s-space 
differential equation 


1 d_ 

s ds 





(40) 


The first term in (l40l) is obtained through integration by 
parts twice, which is only applicable if [64J 

To < -r and ^2 < - as s 0. (41) 

s 2 ds s 

In addition, To(s) must decay sufficiently fast to allow 
convergence of Defining 


Y 0 (s) = e~^ s To(p), p = 2V2s, (42) 


at Z = 0 [we will sometimes write —V for the right-hand- 
side, see (EH)], and the far-field condition 

0 —> 0 as R 2 + 2 2 —y oo. (34) 

The logarithmic dependence of the eigenvalues A and 
their corresponding eigen-potentials 0 upon h is evident 
from (l33l) . We also see that a rough leading-order solu¬ 
tion, entailing an (^(ln -1 L) relative error, may be ob¬ 
tained by neglecting the right hand side of (l33l) . In 
the following subsections we proceed in two paths: (i) 
in 11III El we calculate succesive terms in inverse logarith¬ 
mic powers; (ii) in 11III Fl we tackle the reduced eigenvalue 
problem as a whole, thereby finding an “algebraically ac¬ 
curate” approximation for the resonance e values. 


substitution shows that T(p) satisfies Laguerre’s equation 
d 2 T 0 dTc 


d 2 T° dT 0 1 / a/2 \ , . 

V~TJT + (1 - P>— + o ( — - 1 1 T o = °- ( 43 ) 


•o 


Non-singular solutions of the latter, which are consistent 
with (l4Tj) and the required attenuation rate of To(s), ex¬ 
ist only when the factor multiplying T in (14311 is a non¬ 
negative integer. In those cases the solutions are propor¬ 
tional to the Laguerre polynomials. Hence, the leading- 
order eigenvalues are 


A 


(n) 

0 


V2 

2 n + T 


n = 0,1,2,... 


The corresponding transformed solutions are 


(44) 


E. Two terms in a logarithmic expansion 


T 0 (ti) (s) = /C (n) e-^ s L„(2v / 2s), (45) 


We here expand the potentials 0 an d the eigenvalues 
A as 

0 ~ 0o + t— f 0i + • • • > A ~ Ao + -—-Ai + • • • (35) 

ln h ln h 

At leading order, the boundary condition (l33l) reads 

0o + A o ff(i?)^J =0 at 2 = 0. (36) 

oZ 

We look for a solution in the form 

Ms,z) = n[MR,z)\, (37) 


where 

POO 

H[f(R)](s)= / f(R)J 0 (sR)RdR (38) 
Jo 

is the zeroth-order Hankel transform [64| . Applying the 
Hankel transform to Laplace’s equation (l32l) . a solution 
that satisfies (1M1) is 


00 = -Ao(s)e 

s 


-sZ 


(39) 


where is an arbitrary constant, and L n is the Legen¬ 
dre polynomial of order n. It is straightforward to invert 
(1371) to get the eigen-potentials 0^. For example, not¬ 
ing that Lq(j>) = 1 and L\(p) = 1—p, the first two modes 
are 


4 0) = 


/C(°) 


[R 2 + (Z + v^) 2 ] 


1/2 


and 


(1) (0) 2 v / 2/C«(Z + v / 2) 

Wo = Wo - 


[R 2 + (2 + v^) 2 ] 


3/2 ' 


(46) 


(47) 


Proceeding to the next logarithmic order, we see that 
the problem governing is similar to the one govern- 
ing 0g ', only that the homogeneous boundary condition 
(13611 is replaced by an inhomogeneous one. Dropping the 
( n ) superscript for the moment, the latter condition fol¬ 
lows from (1531) as 


/ m X _ A 1 , ^ f°° ^ R 

0 1 + A o H( J R)^ J --0o + y o 


dR. 


(48) 
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From the preceding analysis it is clear that, when Ao 
is given by m, the problem governing ipi possesses a 
nontrivial homogeneous solution. Hence, the Fredholm 
alternative suggests a solvability condition on the forcing 
terms in (1551) . Indeed, by subtracting (1551) multiplied 
by Ri/jo/H from (1551) multiplied by Ri/ji/H, followed by 
integration with respect to R and use of Green’s second 
identity, we find 


Ai 


—Ao 



H{R) 


dR. 


(49) 


The integrals in (1591) can be readily evaluated by use 
of Parseval’s identity and the inversion-symmetry of the 
Hankel transform, together with the orthonormality of 
the Laguerre polynomials, see ([A.ll) . One finds the sim¬ 
ple result Ai = — 2\/2{X 0 ) , whereby reinstating the (n) 
superscript gives the two-term logarithmic expansion 


A (") 


V2 

2n + 1 



4 1 \ 

2n + 1 ln(l /h) + ) 


(50) 


F. Algebraically accurate result 

It is clear that the logarithmic error in (15U1) limits the 
usefulness of that result to exceedingly small h. To im¬ 
prove on this, we now return to the reduced problem of 
m Dl and attempt a complete solution for the A eigenval¬ 
ues, treating logarithmic terms in par with 0(1) terms. 
This is equivalent to calculating all of the terms in expan¬ 
sion (1551) . We begin by writing the solution as [cf. (1571) . 
(1551) and (1551) ] 


V’ = R 1 [V>], where i/) = -Y(s)e sZ . (51) 

s 

Applying the Hankel transform to (1551) . we find 


l±( dY 

s ds \ ds 



4V Sjs) 
A s 


(52) 


where 6(s) is the Dirac delta function. E q. (l52l) is to 
be interpreted in the sense of distributions [65|, which is 
easier to do when identifying the Hankel transform as a 
two-dimensional Fourier transform of a radially symmet¬ 
ric function, s being the radial distance in Fourier-space. 
Crucially, conditions (1411) are now inapplicable; rather, 
balancing singularities at the origin yields the condition 

2V 

Y(s) ~——Ins as s —> 0. (53) 

A 

Following (1551) . we set Y(s) = e~^ s T{p), where p = 
2\/2 s. This gives [cf. (1551) ] 

cj2j~< dT 

p ~d^ + ^~ p ^ +flT = 0, p> °’ (54 ) 


where h is defined through 


A = 


V2 

2h + 1 


(55) 


We arrive at an eigenvalue problem consisting of (1541) . 
the “boundary” condition 


Tip) 



as p — > 0, 


(56) 


and the condition that T does not grow too fast as 
p —> 00 . The solutions found in ( IIII El which correspond 
to non-negative integer fi and V = 0, are inconsistent 
with condition (1561) . Hence, in the present scheme, where 
logarithmic errors are not tolerated, V / 0. In the lat¬ 
ter case, excluding non-negative integer n, solutions of 
(El) which are logarithmically singular as p —> 0 and 
have acceptable behaviour at large p are proportional 
to U{—n,l,p), where U is the confluent hypergeometric 
function of the second kind [H}. Noting that 

U(-n,l,p) - r ^_. [lnp+ T(-n) + 27 e ] + o(l) (57) 


as p —> 0, where T(a;) is the Gamma function and ^(a;) = 
r , (a;)/r(a;) is the Digamma function, condition (l56l) is 
satisfied by 


T(p) = — r(— n)U(-h,l,p). (58) 

It remains to satisfy the global condition, that is, to 
relate V to the right hand side of (1551) [cf. ([55]) ]. It is 
convenient to first rewrite that condition as 

2 f 00 / (9 1 A 

b s +2 W„ (59 > 


The integral on the right can be evaluated by using Parse- 
val’s identity. Noting that the Hankel transform of 1 /H 
is 2A 0 (sv / 2), K q being the modified Bessel function of 
the second kind, that integral is 


Given ([57l) . and since K 0 {x) ~ — ln(x/2) — y e as x —> 
0, the apparent logarithmic singularity of the brack¬ 
eted terms cancels, whereby the sifting property of the 
Dirac delta function can be effected to show that © 
= 2 In4 + 2'F(— n) + 2 r ) e . Substitution of this result into 
(1591) furnishes a transcendental equation for h: 



-h-Y(s) + 2K 0 (sV2 


6{s)ds. (60) 


2'F(—n)=lnA. (61) 

Since ^(z) diverges at non-positive integer x, it is easy to 
see that ED has multiple non-integer solutions fi which 
(very slowly) approach n = 0,1, 2, as h —>0 [67j] . 


















FIG. 4. The radial distribution of the dimensionless trans¬ 
verse field in the gap for the longitudinal modes n = 0,1 and 
2 , for h = 0.001 and V = 1 . 


Recapitulating, definition (1551) furnishes the eigenval¬ 
ues as 


a(") 


V2 

2 h(h) + 1 ’ 


(62) 


where h is the solution of (ED which approaches n = 
0,l,2,...ash —► 0. We reiterate that (l62l) is algebraically 
accurate, viz., together with m it provides the reso¬ 
nance values of the relative dielectric function up to a 
relative error which scales with some power of h. In par¬ 
ticular, (1621) is far more accurate then (l50l) . We note that 
it is straightforward to recover (1501) by expanding tit in 
ED about its poles. Higher order logarithmic terms can 
be derived in a similar manner. Inverting if) provides the 
distribution of the eigenpotentials in the pole to algebraic 
order [cf. (HtH) and (ITfl) ]. The dominantly transverse field 
in the gap, which is a function of R alone, is then obtained 
(up to a minus sign) as d<&/dZ = A(R) = 0)/H(R). 

The radial gap-field distributions of the first three modes 
are plotted in Fig. (SH¬ 


IV. INDUCED-DIPOLE MOMENT AND FIELD 
ENHANCEMENT 

A. Resonance conditions 

We now return to the excitation problem formulated 
in where the far-field condition © applies, and e is 
a complex-valued function of frequency. Whereas a neg¬ 
ative real e is unphysical, the reduced Drude model © 
implies that — Re[e] Im[e] as long as w > 7 . Reso¬ 
nance frequencies are safe within this range: even for an 
extremely narrow gap, say h = 10 -4 , Eqs. © and CD 
suggest uj ~ 0.1 w p — still quite larger than 7 for silver 
and gold. Thus, towards studying the damped near-field 
resonance of the system, it is useful to consider a reso¬ 
nance scenario 

e ~ e(") + iti, €i <C (63) 

Here, for a given n = 0,1, 2,..., e( n ) is a resonance value 
~ —A( n )(/i)/i -1 / 2 + 0 ( 1 ), where A( n ) = 0 ( 1 ) is given by 


(1621) . and iti is the imaginary deviation from the reso¬ 
nance values. 

Given (1631) . it turns out that an asymptotic descrip¬ 
tion of the field distribution — and in particular, closed- 
form approximations for the induced-dipole moment 7 x 
and field enhancement in the gap G — can be obtained 
based on the LSP modes calculated in m The crucial 
observation is that the nth LSP mode is “excited” un¬ 
der (1631) . viz., the leading-order potential is of the form 
calculated in m but with an asymptotically large mul¬ 
tiplicative factor. Note that if the potential was not 3> 1, 
(RJ51) would imply a contradictory leading-order problem 
which, on the one hand, possesses a non-trivial homoge¬ 
neous solution, and on the other, is forced at large dis¬ 
tances by the incident field. Thus, the integral relation 
© can be written, asymptotically, as 

47/ V 79 ■ V 93 d.V. (64) 

J metal 

We can now substitute our solutions from ( 1IIII (allowing 
for a large pre-factor) into (l64l) . and solve for the am¬ 
plification factor (say, in terms of /x). We note that, in 
contrast to ED: the work-dissipation relation © is am¬ 
biguous with respect to phase and hence cannot be em¬ 
ployed towards making quantitative assessments without 
further information. We once more precede a detailed 
analysis with scaling arguments. 


B. Scaling arguments 


Let <fg denote the order-of-magnitude of the poten¬ 
tial in the gap region. The asymptotic structure of the 
potential then follows from the description provided in 
m In particular, the transverse field in the gap is 
G = 0(ip g /h), the potential in the pole region is 0(ip g ), 
and on the particle scale the sphere is approximately 
equipotential, with V = 0((p g / ln(l//i)), see (13T1) . Note 
that ED implies /x = O(V). Consider now the integral 
on the right hand side of ( 1 M 1 ) [the work-dissipation re¬ 
lation © can also be employed]. The dominant contri¬ 
bution arises from the pole domains, where the potential 
varies rapidly. Indeed, the volume of the pole domains is 
0(h 3 / 2 ) and the integrand there is 0((p 2 g /h ), leading to an 
O^g/i 1 / 2 ) estimate. The contribution from the approx¬ 
imately equipotential 0 ( 1 ) volume is 0(ip 2 h) at most. 
We thus find from (l64ll that ip g = 0(/i - 1 / 2 /(e* ln(l//i))), 
from which estimates for /x and G readily follow as 


O 


1 


: h 1 / 2 In 


2 1 


G~0 


G h 3 / 2 In i 


(65) 
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C. Asymptotic formulas 

To go beyond scaling arguments, we use Green’s first 
identity together with © to rewrite (1551) as 


To obtain algebraically accurate formulas, we substi¬ 
tute (1551) for Y and (1551) for A. Making use of the in¬ 
tegrals (1A.3I) and (1A.4I) . the solution scheme described 
above yields 


fj, ~ -ieih 1/2 J _ RdR i (66) 

where denotes the deviation of the potential in the 
(z > 0) pole domain from V, see dllll Applying Parse- 
val’s theorem to the right-hand-side of (1551) . it becomes 

[cf. m] 


pOO 

ie x h 1/2 / Y 2 (s)ds. 

Jo 


(67) 


Upon substitution of Y (see below), and noting that \i ~ 
7t 2 V/3 [cf. (1571) ]. (1551) can be solved for V, and hence /i. 
The field enhancement we can then be obtained from 


G 


1 d<f» 
h~dZ ’ 


( 68 ) 


V^TT 4 _ 

9iej/i 1 / 2 (2n + l) 2 d>'(—n) 


and 


v^tt 2 T(—n) 2 Fi (2,2,2 - n, 1/2) 
12ieih 3 / 2 (2h+l) 2 4"(-n)T(2-h) ’ 


(73) 


(74) 


where ^'(j;) is the TriGamma function, the derivative 
of the DiGamma function, 2 Fi is the hypergeometric 
function, and h is as before the solution of © that 
approaches n = 0,1,2,... as h -A 0. We note that 
(ED and ED can be derived as leading approximations 
of ed and ED. respectively, by noting that n — n ~ 
2/ln(l /h) as h — > 0, and that 4 ,, (— fi) ~ — (n — n) -2 , 
2 Fi(2,2,2-n,l/2)/r(2-n) ~ 4n!(2n+l) and T(-n) ~ 
(— l)" +1 /[n!(n — n)] as n —> n. 


where $ is the gap potential [cf. ED], and the derivative 

is evaluated at R,Z = 0. Using ED, G can also be V. COMPARISON WITH SEMI-NUMERICAL 
written as SOLUTION 


\dj}_ 

h~dZ' 


(69) 


where the derivative is evaluated at R = 0, Z = 0, or in 
terms of Y (s), 


A r°° 

~-t/ sY(s)ds. (70) 

hj o 

The calculation can be performed at two levels, de¬ 
pending on whether the logarithmically or algebraically 
accurate expressions for Y and A are employed in the 
above formulas. We start with the former case, substi¬ 
tuting 1q(s) f° r Y ( see flip ) and Ao for A (see (1551) '). 
Noting ED, the arbitrary constant K, in Vq is to within a 
logarithmic relative error /C ~ V In //Ao. Together with 
the integrals (IA.1I) and (I A. 2D , we find the explicit closed- 
form expressions 


4\/27r 4 

^ 9(2n + l) 2 iei h 1 / 2 ln 2 (l/ft.) 

and 


(71) 


(-1 ) n + 1 2 v / 27T 2 
3(2n + 1 )iti h?! 2 ln(l/ h) 


(72) 


Expressions ED and d72l) have the advantage of being 
simple, but since they entail an 0(1/ ln(l//i)) error, they 
are hardly accurate, as we shall later demonstrate. We 
note the agreement with the scaling results (1551) , and that 
the excited resonant mode is imaginary, i.e. out of phase 
with the incident field. 


We have validated our asymptotic formulas against a 
semi-numerical solution of the quasi-static problem. We 
employ a scheme which is based on the separability of 
Laplace’s equation in bi-spherical coordinates. Briefly, it 
consists of writing down a formal solution in the form 
of two infinite series, and determining the coefficients by 
numerically solving a truncated tri-diagonal system of 
coupled algebraic equations. The relevant formulas are 
available elsewhere (*'>. H, H(| . 

An example calculation is shown in Fig. [D where we 
have employed the empirical e-data for silver provided by 



1.5 2 2.5 3 3.5 

eV 


FIG. 5 . Absolute magnitude of the dimensionless induced- 
dipole moment fj, as a function of energy hu and h, half the 
dimensionless gap width. The coloured contours represent 
the semi-numerical solution of the quasi-static problem. The 
enumerated white lines depict the resonance frequencies of 
the first three modes as predicted by the asymptotic formula 
( 1621 ) . together with ED- We employ the empirical silver data 
of Johnson and Christy [ 68 :]. 
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FIG. 6. Zeroth-mode e-eigenvalue as a function of h, half the 
dimensionless gap width. Results of a semi-numerical solution 
(symbols) are compared with the first term and two terms of 
the logarithmic expansion m, and the algebraically accurate 
approximation (1621) . cf. (1101) . Note also the approach of e to 
—2 at large h. 

Johnson and Cristy 0]. Fig. [5] shows the absolute mag¬ 
nitude of the induced-dipole moment g as a function of 
the dimensional energy huj in eV units, for a range of h. 
We observe magnitude peaks which, as expected, redshift 
as h —> 0. The white lines in Fig. [5] depict the first three 
resonance frequencies as predicted from the asymptotic 
formula (1621) in conjunction with the data of Johnson and 
Cristy. For small h the latter lines successfully trace the 
amplitude peaks. Detailed results for the zeroth-mode e- 
eigenvalues are shown in Fig. [6] It depicts as a function of 
h the numerically obtained eigenvalues (symbols) along 
with the predictions of (1621) and the two-term expansion 
(l50l) . The asymptoticness of (l62l) as h —> 0 is evident, 
with good agreement for h < 0.02 (i.e. radius-to-gap ra¬ 
tio of ~ 25). As expected, the logarithmically accurate 
approximations are quite poor. Note the approach at 
large h of e to — 2, the Frohlich value of isolated spheres. 

Consider next the values of /i and G. Focusing again on 
the zeroth-mode, we have recorded for a range of h the 
numerical values attained by fx and G at their zeroth- 
mode peaks. The absolute magnitude of these values 
(which are predominantly imaginary at resonance) are 
depicted by the symbols in Fig. [3 The latter numeri¬ 
cal values are compared with the logarithmically accu¬ 
rate asymptotic approximations ED and ED — dashed 
lines, and the algebraically accurate approximations ED 
and ED thick lines. Again we see that the latter 
approximations are rather good, whereas the former are 
poor for realistically small h. 

VI. DISCUSSION 

We have employed singular perturbation methods to¬ 
wards analyzing the longitudinal surface modes of a pair 
of metal spheres in near contact, and the quasi-static 
excitation of these by an electromagnetic plane wave. 
Whereas the techniques of matched asymptotic expan¬ 



FIG. 7. Numerical peak values of the dimensionless induced- 
dipole moment fi and the field enhancement in the gap G at 
zeroth-nrode resonance, as a function of h, half the dimension¬ 
less gap width. Also shown are the logarithmically accurate 
approximations ED and ED. and the algebraically accurate 
approximations ED and ED- We employ the empirical silver 
data of Johnson and Christy [68]. 

sions are routine in many areas of research, it appears 
that their application in the present context is novel. 
As demonstrated herein, these allow describing limits of 
physical extremity in a relatively simple and intuitive 
manner, notably without the prerequisite of an analyti¬ 
cally cumbersome or computationally heavy “exact” so¬ 
lution. In particular, the present analysis has furnished 
accurate asymptotic expressions for the resonance values 
of the dielectric function (or frequencies) — see (1621) and 
ED- and the values at resonance of the induced dipole 
moment and field-enhancement in the gap - — see ED and 
(T74l) . These expressions have been validated in ifVlagainst 
a semi-numerical solution of the quasi-static problem. 

The above-mentioned expressions are asymptotic in 
the limit h —> 0, with algebraic accuracy. That is, they 
involve a relative error on the order of some power of 
h. These results were preceded with simpler expres¬ 
sions — see ED. ED and ED. which incur however 
a much larger relative error on the order of l/ln(l/h). 
As the comparison in [JV] unsurprisingly suggests, the lat¬ 
ter logarithmic expressions are far less accurate for rea¬ 
sonably small h. Nevertheless, the analysis leading to 
these revealed the essential structure of the longitudi¬ 
nal LSP modes, thereby laying a path towards their al¬ 
gebraic counterparts. We note that similar logarithmic 
expressions have been derived by approximate manipu¬ 
lations of implicit series solutions in bi-spherical coordi¬ 
nates [H, j5(|. The limited value of these expressions, 
however, was not emphasised in those works. To be spe¬ 
cific, the leading term in (l5(JD agrees with both references, 
and (1721) agrees with the expression given in [50:]. The 
logarithmic correction for e given in Ref. |o(| , however, 
appears to posses a factor-2 error when compared with 
the second term in (l5UP . 

In addition to furnishing closed-form approximations, 
our analysis offers insight into the physical and geometric 
circumstances enabling a strong frequency redshift and 
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field enhancement. The former has to do largely with 
the gap morphology, which enables stronger field local¬ 
ization in the dielectric than in the metal — a prerequi¬ 
site for large-1e| eigenvalues. In fact, the scaling of e with 
1/h 1 / 2 was shown in £1111 Al to follow from the particular 
imbalance in transverse localization pertinent to a locally 
parabolic boundary. This scaling should therefore apply 
quite generally to the large-|e| resonances of essentially 
arbitrary plasmonic two- and three- dimensional clusters 
in near contact, or to particles in near contact with a 
substrate (see Ref. [13 for an alternative intuitive argu¬ 
ment). The origin of field enhancement is more subtle. 
The amplitude of an excited mode, represented e.g. by 
the configuration-scale induced-dipole moment, is deter¬ 
mined by a balance of work and dissipation, which is in 
turn affected by the geometrically enabled localization. 
Disregarding a weak logarithmic dependence upon h 7 the 
latter enhancement is on the order of |e/e;|, which in the 
near-contact limit is 0(/i _1 / 2 /ei). The giant 0(h~ 3 ^ 2 /ei) 
enhancement in the gap represents a further 0(l/h) geo¬ 
metric amplification relative to the induced particle-scale 
field. We note that in the corresponding planar problem 
of a nanowire dimer, where scaling arguments similar to 
those employed herein reveal that the longitudinal modes 
are determined solely from the gap morphology, the en¬ 
hancement is smaller, 0(h~ 1 /ei) [43]. This fundamental 
difference is apparently overlooked in Ref. [§|, where in the 
case of a spherical pair the resonant field enhancement 
in the gap is estimated incorrectly on pg. 278 [cf. (l65jl ]. 

For the sake of demonstration, we have focused on the 
axisymmetric longitudinal modes of the spherical dimer 
configuration, which strongly redshift with vanishing sep¬ 
arations and are the ones excited by a field polarized 
along the line of sphere centers. In many applications 
the latter modes play the key role, but this is not al¬ 
ways the case, and it is therefore desirable to apply the 
present approach to LSP modes having different sym¬ 
metries (48(. For example, it has been shown that non- 
axisymmetric modes of high azimuthal order provide the 
dominant contribution to van der Waals (vdW) energies 
at small separations [70]. In Ref. [?0, approximate expres¬ 
sions for the vdW energies of sphere dimers and a sphere 
near a wall were derived by reduction of implicit infinite- 
series solutions in the double limit of small separations 
and high azimuthal mode number. While, as already dis¬ 
cussed, similar methods gave in Ref. 51 only logarithmic 
accuracy in the case of longitudinal modes, the dominant 
modes in the vdW calculation are more highly localized 
in the gap and accordingly do not involve a logarithmic 
coupling with the particle-scale field. 

Our methodology can be applied to study various other 
plasmonic nano-structures that are characterised by mul¬ 
tiple length scales. In particular we note the related 
particle-substrate configuration, which can be addressed 
along the lines of the present analysis, the substantially 


simpler two-dimensional problem of a nanowire dimer, 
and elongated particles. Further desirable generalisa¬ 
tions are concerned with physical modelling. Our anal¬ 
ysis originates from a quasi-static formulation, and fo¬ 
cuses on the near-field response of a small nano-metric 
particle (still large compared with the gap width). Given 
our expressions for the induced-dipole moment, estimates 
for the optical cross section can be obtained via the 
usual connection formulas @, Qjj, HU . Ad hoc techniques 
for improving and extending such estimates to larger 
nano-structures, wherein the quasi-static induced-dipole 
is renormalised to account for radiation damping, have 
proven quite useful jUIzJ ■ More rigorous however would 
be to connect the near and far fields through systematic 
application of matched asymptotic expansions [65l . [72j , 
taking special account of the largeness of |e|, the large¬ 
ness of the near-field at resonance, and the multiple-scale 
geometry of the metallic structure; strictly speaking, each 
of these disparities entails slaving a distinct small pertur¬ 
bation parameter to the Helmholtz parameter aw/c. 

Probably the most important direction is to incorpo¬ 
rate into the present study non-classical physics, which at 
sub-nanometric separations ultimately su ppr ess the fre¬ 
quency redshift and field enhancement [29j, l45l . l73| _. In the 
framework of the hydrodynamic Drude model [74L [75}, 
my co-workers and I recently showed [76 ] that the near¬ 
contact asymptotics predicted by the latter “nonlocal” 
model constitute a renormalisation of those predicted by 
the classical “local” model considered herein. This find¬ 
ing greatly extends the utility of the present study. 
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Appendix: Integrals 

We employ the integrals 

pOO 

/ e~ p [L n (p)] 2 dp = 1, (A.l) 

Jo 

poo 

/ pe~ p ^ 2 L n (p) dp = 4(—1)"(1 + 2n), (A.2) 

J o 


1 o T'(-fi) 

e- p [U(-h 7 l,p)} 2 dp= 1 j 


(r(-n)r 


(A.3) 


pe p l 2 U{—h, l,p) dp = 


2^1 ( 2 , 2,2 — h , 1 / 2 ) 


F(2 - n) 


(A.4) 


for n = 0,1,2,... and non-integer n > 0 
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